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Interactions of dispersive shock (DSWs) and rarefaction waves (RWs) associated with the Korteweg-de Vries 
equation are shown to exhibit multiphase dynamics and isolated solitons. There are six canonical cases: one is 
the interaction of two DSWs which exhibit a transient two-phase solution, but evolve to a single phase DSW for 
large time; two tend to a DSW with either a small amplitude wave train or a finite number of solitons, which 
can be determined analytically; two tend to a RW with either a small wave train or a finite number of solitons; 
finally, one tends to a pure RW. 
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Shock waves in processes dominated by weak dispersion 
and nonlinearity have been experimentally observed in plas- 
mas |U), water waves [2|, and more recently in Bose-Einstein 
condensates (3] |4) and nonlinear optics [5|; these dispersive 
shock waves (DSWs) have yielded novel dynamics and inter- 
esting interaction behavior which has only recently begun to 
be studied theoretically (cf. EH). Here we consider DSWs 
which are described by the Korteweg-de Vries (KdV) equa- 
tion, 



u t + uu x + s u xxx = 0, < s <k 1 . 



(1) 



Individual DSWs are characterized by a soliton train front 
with an expanding oscillatory wave at its trailing edge; these 
waves have been well-studied (cf. |8, 9]) using wave averag- 
ing techniques, often referred to as Whitham theory IfTOlfTTI . 

When illustrative, we contrast DSW interaction with classi- 
cal or viscous shock waves (VSWs), which are dominated by 
weak dissipation and nonlinearity, using Burgers' equation 



u, + uu x — vu xx = 0, 0<v« 1. 



(2) 



The interaction of VSWs is an entire field and has been ex- 
tensively studied (cf. [ 12 1), while little is known about DSW 
interactions. 

In this letter, we use analytic, asymptotic and numeric 
methods to investigate ([T]) (and using the "step-like" ini- 
tial data 



u(x, 0) = Uq(x) 



ho, x < 0, 
h\, < x < L, 
hi, x > L, 



(3) 



where ho, hi and hi are distinct, real and non-negative. This 
gives six canonical cases, which we denote: 




ho > hi > hi, 
hi > ho > hi, 
hi > hi > ho, 




): ho> h 2 > h u 
): h 2 > ho > hi, 
): h 2 > hi > ho, 



where an icon of the initial step data is shown in parentheses. 
When convenient, and without loss of generality, we take /i, 
to be 0, 1 and < h, < 1 (by using a scaling symmetry and 
Galilean invariance). The case of an initial depression (e.g. 
Case II, ho = hi = > hi) and an initial box (e.g., Case 
III, ho — hi — < hi) has been studied in [7|, where the 
asymptotic solution was constructed analytically. 

This letter is organized as follows. We first discuss Case I 
C~ l -lJ, where two DSWs interact and exhibit a two-phase re- 
gion which evolves into effectively a one-phase solution for 
large time. Single phase Whitham theory is then introduced 
to describe the DSW with a small amplitude wave train which 
develops in Case II (H_i— )■ We then briefly discuss multi- 
phase Whitham theory to describe the two-phase region in 
Case I (~~ •— u_)- In Case III (— r l_), the interaction produces a 
DSW with a finite number of solitons, which remarkably can 
be determined analytically using Inverse Scattering Transform 
(1ST) theory (cf. |[T8l ). There is no analogue for emerging 
solitons in VSWs. We then use Whitham and 1ST theory to 
describe the interactions in Case IV (— i_p ), V (_P— ) and VI 
(_ I — r ~ ). Finally, we comment on the numerical scheme we 
used to solve ([TJ and |2]). 

In Case I ( — 1 ^_), two one-phase DSWs form and propa- 
gate to the right (see Fig. When the shock front of the 
left DSW reaches the expanding oscillatory tail of the right 
DSW, they interact and form a quasi-periodic two-phase so- 
lution (see Fig. [TJi). The shock front of the left DSW sub- 
sequently overtakes the shock front of the right DSW and 
forms a one-phase solution to the right of the two-phase re- 
gion (see Fig.[T]:). To the left of the two-phase solution, an es- 
sentially one-phase DSW tail emerges (see Fig.[T];); although 
the tail is weakly modulated by a quasi-periodic wave, its be- 
havior is essentially one-phase. For large time, the two-phase 
region closes and a one-phase DSW remains (see Fig. [TJl- 
e); Whitham theory indicates that the amplitude of the two- 
phase modulations decrease with time and result in an effec- 
tively one-phase DSW. This closing of the two phase region 
is suggested by the rigorous (Whitham theory) results in lfT4l . 
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FIG. 1: Plots (a)-(d) show the numerically computed solution of Q 
and (e) the boundary of the one- (light gray) and two-phase (dark 
gray) regions computed using Whitham theory. The averaged solu- 
tion, u, is computed using Whitham averaging (cf. 1 13 1) and shown 
as dotted lines in (a)-(d); the solution of (|2j is shown as dashed lines 
in (a)-(d). In all plots, s 2 = 0.001, h = 1, hi = 0.4, h 2 = and 
L = 8. The vertical axis in (e) is log-time and the horizontal axis is 
—t < x < t + 8 (and matches the domain in (a)-(d)). 
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FIG. 2: Plots of the numerical and averaged Whitham solutions of 
(]} for Case II n_r-) where b 2 = 0.001, ho = 1, hi = 0, h 2 = 0.5 
and L = 5. 



though the authors studied smooth initial data. The computa- 
tion of the boundaries of the one- and two-phase regions using 
multiphase Whitham theory are discussed later in this letter. 

Although the (initial) shock front speed is different for 
DSWs and VSWs (2/io/3 and ho/2, respectively), the aver- 
aged DSWs are similar in behavior to VSWs (see Fig.[T^-d); 
in both, two shock waves merge to form a single shock wave. 

For Case II C~l__r- ), a large DSW forms on the left and a 
small RW forms on the right (see Fig. |2^). The front of the 
DSW then interacts with the trailing edge of the RW; the in- 
teraction decreases the DSW's speed and height (see Fig. [2}}). 
The front of the DSW is faster than the front of the RW and 
overtakes it (see Fig. |2f;). The size of the interaction region 
continues to expand with a DSW emerging in front with a 
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FIG. 3: The initial data regularization of Case II ( | , — ) for h > 1, 

hi = and hi = 1; the dashed line is the initial condition, u (x), and 
the solid lines are r\, r 2 , and 73. The figure also gives the speed of 
the front and back of the DSW and RW at / = 0. 



small amplitude wave train behind, whose amplitude is pro- 
portional to r 1 ^ 2 (see Fig.[2jl). As in Case I ( — ,_), the av- 
eraged DSW and the VSW (see Fig. |2]i both tend to a single 
DSW (VSW) once the front of the DSW (VSW) passes the 
front of the RW. 

We can use the one-phase Whitham equations to character- 
ize the interaction of the DSW and RW in Case II ( — ). In 
this context, Whitham theory consists of looking for a fully 
nonlinear single- or multi-phase solution whose parameters 
(amplitude, wave number and frequency) are slowing varying 
with respect to the phase(s) and then deriving new equations 
for the evolution of the slowly varying wave properties. The 
one-phase Whitham equations for ([l} are 



— + Vi(n,r 2 ,r 3 )— = 0, 
at ox 



i = 1,2,3, 



(4a) 



where 



V- 


3 (r2 


-ri) 


V- 


~(r 2 


-n) 


v + 




-n) 



K(m) 



K(m)-E(m) 
(1 - m)K(m) 

E(m) - (1 - m)K(m) ' 

(1 -m)K(ni) 
E(m) 



(4b) 



V = (ri+r2+7"3)/3, m = (r2-ri)/(r^-r\), K(m) is the complete 
elliptic integral of the first kind, and E(m) is the complete el- 
liptic integral of the second kind (SJ. Then, the asymptotic 
solution is 

u a (x, f) « r\ + r% — r-i + 2(ri - r\) dn 2 (6; m) , 



where Q x = k, 6 t = -to = —kV, k = V( r 3 _ r i)/(6e 2 ), and r, 
are slowly varying functions of x and t. We can make a global 
dispersive regularization for the initial value problem ([TJ and 
Q by choosing appropriate initial data for the r, ||3][15] which 
result in a global solution. A global dispersive regularization 
of Case II (~ l_r-) is shown in Fig. |3j the r, are taken to be 
nondecreasing, r,(x, 0) < r, + i(x, 0) and u„(x, 0) = u(x, 0) for 
all xeR. 

In order to study the interaction we evolve the r,- numeri- 
cally. A simple and effective method for evolving the r, is to 
discretize the initial data regularization along the dependent 
variable, r,, and then compute the shift in x of each data point 
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FIG. 4: Plot (a) shows the Whitham approximation and (b) direct 

numerics of the solution of (|TJ for Case II ( | , — ) with the same 

initial condition as Fig. [2] 



using Q. Fig. |4] compares a numerically evolved Whitham 
approximation with direct numerics for Case II (~ l_r— ); the 
first order Whitham approximation does not capture the small 
quasi-periodic modulations in the tail because they are higher 
order effects. Both direct numerics and the Whitham approxi- 
mation agree and show that for large enough time, the ampli- 
tude of the tail in Cases II ( — ) is proportional to r 1 ^ 2 ; this 
is typical of a uniform linear wave train when the total energy 
remains constant (cf. 1 10 1) and was observed in the context of 
a depression initial condition in |7 1. 

Multiphase Whitham theory is more complicated than one- 
phase Whitham theory and dates back to 1970 [16|; multi- 
phase Whitham equations were developed for the KdV equa- 
tion in IfTTl . The interaction of two DSWs from certain 
step-like data was recently analyzed in [ 6 1 for the nonlinear 
Schrodinger equation. The one- and two-phase regions and 
the averaged solution in Case I (~~ •— ,__) are found by numeri- 
cally evolving the two-phase Whitham equations for the KdV 
(see US), 



On 
dt 



,r 5 )^ = 0, 
ox 



i = 1,2,..., 5, (5) 



where v,- = (2r] ~xr 2 -Bm-fS^/ir 2 - am- a 2 ),X = E,=i r j 
and a\, a 2 , B\ and B 2 are solutions of 



j\ jb 

L 2 L 2 



1 1 l 
1 2 1 2 



with 



(6) 



In Case III (— H_), a small RW forms on the left and a large 
DSW forms on the right. The front of the RW then interacts 
with the tail of the DSW and reduces the amplitude of the 
waves — essentially cutting off the top of the box. Since the 
front speed of the RW is less than the front speed of the initial 
DSW, a finite number of solitons can escape the interaction 
(see Fig. [5]). These solitons have no analogue in the VSW 
solution of Case III (— r ~l^. We can compute the number of 
solitons which escape using 1ST theory. 

From 1ST theory, the number of solitons correspond to the 
time-independent number of zeroes of a(k) (which is the num- 
ber of poles of the reflection coefficient R = b(k)/a(k)) in the 
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FIG. 5: Plots of Cases (a) & (b) III (-H_) with h = 0.5, hi = 1, 
h 2 = and (c) & (d) V (_P— ) with h = 0, h = 1, h 2 = 0.5, where 

b 2 



upper half £-plane. Associated with ([TJ, the data a(k) is de- 
fined by 

(f>(x\ k) = a{k)\p(x; k) + b{k)\]/{x; k), 
0(jc; k) = a{k)\p(x; k) + b{k)\j/{x; k), 

corresponding to the eigenfunctions, 



4>(x; k) ~ e- ikoX , <f>(x; k) 

i//(x;k) ~ e' klX , i//(x;k) ~ e mz % as x — > -fco. 
which satisfy the Schrodinger scattering problem, 

0. 

The solution of (|7), at t = 0, is 

' Ae ik x + B e - ikoX , x < 0, 
w(x) = | Ce ihx + De~ iklX , < x < L, 
Ee ik * x + Fe~ iklX , x> L, 



e ikoX , 

-ikix 



as x ■ 

as x 



W xx + w{u/6 + k 2 }/s 2 



(7) 



where k u = + k 2 /s, k\ = + k 2 /s, and k 2 = 

/s. The eigenfunctions, (p, cp, ip- and \p are deter- 
mined by requiring that w and w' are continuous across x — 
and x = L. Indeed, <p is found by taking A — and B — 1 and 
then solving for C, D, E = b(k), F = a(k), so that 



2k 2 I k l (k + k 2 ) ) 

Since e' klL (k u + k 2 )/(2k2) + 0, the zeroes of a(k) occur when 
tan(feiL) = iki(k + k 2 )/(k 2 + k k 2 ). 

It can be shown that the zeroes of a(k) are purely imaginary; 
thus, we let k = k (where k € R and k > 0). For Case III 
(— r l_J, where h\ - 1 > h v = h* and h 2 = 0, the zeroes of 
a(ii<) occur when 



, , , ^\I6-k 2 (^k 2 -KI6 + k) 

tan ( y/l/6-K 2 L/s) = > ' . (8) 

v ' 1/6 -K 2 -K^K z -hJ6 

The number of periods for y/h*/6 < k < y/l/6 of the RHS 
of flSb, L VI - fh/(sn V6), is an estimate of the number of 
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FIG. 6: Plots of the solution of {J} for Case IV ( — i | ) where e 2 = 

0.001, h = 0.5, h l = 0, h 2 = 1 and L = 5. 

solitons. The number of zeroes determined using (|8]l exactly 
corresponds to the number of solitons observed using direct 
numerics (for various values of /z„, L and s)l 

In Case IV (— i_P ), a small DSW forms on the left and 
a large RW forms on the right (see Fig. |6^). As in Case II 
C~\_!—), the front of the DSW interacts with the trailing edge 
of the RW and decreases the DS W's amplitude and speed. Un- 
like Case II (~~ ), the front of the DSW does not overtake 
the front of the RW. The DSW becomes a small amplitude tail 
on the left of the RW and decreases in amplitude proportional 
to r 1/2 (see Fig. [(§>). 

As in Case III (— H_), Case V (_P— ) cannot be completely 
characterized using Whitham averaging. For Case V (_P— ), 
a large RW forms on the left and a small DSW forms on the 
right; the front of the RW interacts with the tail of the DSW 
and results in a RW and a finite number of solitons. The num- 
ber of solitons corresponds to the number of zeroes of 
where where ho = and hi = 1 > hi — h*. 

In Case VI (_ ! — r ~ ), two rarefaction waves form; the small 
amplitude oscillatory tail (see for instance the RW in |6^) of 
the right RW interacts with the front of left RW; the tail of 
the right and left RW then interact to form a small amplitude, 
modulated, quasi-periodic tail; this modulation decreases with 
time and Case VI (_ I — r ~ ) tends to a pure RW for large time. 

We numerically solve ([T]i and Q using an adaptation of the 
modified exponential time-differencing fourth-order Runge- 
Kutta (ETDRK4) method (see flU]). We use this (sophisti- 
cated) numerical method because ([T]) is very stiff and standard 
numerical methods require the time step to be 0(s 3 ), while for 
ETDRK4 the time step need only be (9(e). When this numer- 
ical scheme was used to compute a known exact solution, it 
was accurate to more than six decimal digits. 

For spectral accuracy when using the ETDRK4 method, the 
initial data must be both smooth and periodic. Therefore, 
we differentiate ([T| with respect to x and define v = u x to 
get v t + (uv) x + erv xxx - 0. Transforming to Fourier space 
gives v, = is 2 k 3 v - ikuv = Lv + N(v, f), where we de- 
fine (Lv)(fc) = is 2 k 3 v and N(v,f) = N(v) = -HcT{[h + 
£ m T~ x (v)dx'YF~ x (v)}. It is important that the integral in N 
is computed using a spectrally accurate method. Moreover, 
we approximate the initial step data with the analytic function 
2wv(x, 0) = (h 2 - hi) sech 2 [(;c - L)/w] + (hi - /i ) sech 2 (x/vv), 
where w is small. See flT9j for details about how this L and N 
are used to numerically compute the solution of ([T]). 



For large time Case I (~ 1_) and II (~|_r— ) go to a sin- 
gle DSW, while Case IV (— uJ~) and VI (_ s - r l go to a sin- 
gle RW; this is consistent with VSW theory. However, un- 
like VSW theory, Case III (-H_) and V (_P— ) form a finite 
number of solitons in addition to the DSW or RW, respec- 
tively. Moreover, unlike VSW theory, Case I ( 1 — i_J exhibits a 
transient two-phase region and Case II (~~L_i— ) and IV (— U - ) 
have a small amplitude tail which decays at a rate proportional 
tor'/ 2 . 
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